package scu.maqiang.cte;

import scu.maqiang.mesh.ScalarFunc;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.fes.FES2Q42;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class RectangleResultQ4_SC1 {

	public static void main(String[] args) {

	    double kxx = 1.0, density = 1.0, specificHeat = 1.0;
	    double E = 2.0, Nu = 0.3, alpha = 0.05;
	    double beta = E / (1 - 2 * Nu)  * alpha;

	    double theta = 0.5, endTime = 1.0e1, dt = endTime / 1000;

	    int verbosity = 1;
	    
	    String XYStressYFile;

	    Mesh2Q4 mesh2D = new Mesh2Q4().square2D(10, 20);
	    mesh2D.scale(0.5, 1.0);

	    //定义热传导问题有限元空间以及相关刚度矩阵与质量矩阵
	    FES2Q41 fs = new FES2Q41(mesh2D);
	    SRMatrix A = new SRMatrix(fs.GetNdof());
	    SRMatrix M = new SRMatrix(fs.GetNdof());
	    double[] x = new double[fs.GetNdof()];
	    double[] RHS = new double[fs.GetNdof()];
	    
	    fs.assembleStiff(new double[] {kxx}, BVPType.COMMON, A);
	    fs.assembleMass(new double[] {density * specificHeat}, BVPType.CONSIST_MASS, M);
	    
	    SRMatrix LMatrix = M.add(theta * dt, A);
		SRMatrix RMatrix = M.add((theta - 1) * dt, A); 
	    
		FES2Q42 fs2 = new FES2Q42(mesh2D);
		SRMatrix K = new SRMatrix(fs2.GetNdof());
		double[] d = new double[fs2.GetNdof()];
		double[] heatLoad = new double[fs2.GetNdof()];
		double[] totalLoad = new double[fs2.GetNdof()];
		double[][] disp = new double[2][fs.GetNdof()];
		double[] ex = new double[fs.GetNdof()];
		double[] ey = new double[fs.GetNdof()];
		double[] exy = new double[fs.GetNdof()];
		double[] tx = new double[fs.GetNdof()];
		double[] ty = new double[fs.GetNdof()];
		double[] txy = new double[fs.GetNdof()];
		
		fs2.assembleStiff(new double[] {E, Nu}, BVPType.PLANESTRAIN, K);
	    
	    fs.applyBC_MBN(LMatrix, Direct.All, 3);
	    fs2.applyBC_MBN(K, Direct.X, 2, 4);
	    fs2.applyBC_MBN(K, Direct.Y, 1);
	    
	    NewIterSSolver solver = new NewIterSSolver(LMatrix);
	    NewIterSSolver solver2 = new NewIterSSolver(K);
	    ScalarFunc betaFunc = (xy, label, param) -> beta;
	    ScalarFunc[] heatLoadFunc = {betaFunc};
	    double[][] vecX = new double[1][];
	    
	    String dir = "RectangleResultQ4\\";
	    int lineNum = 100;
	    double dis = 1.0 / lineNum;
	    double[] lineX = new double[lineNum + 1];
	    double[] lineY = lineX.clone();
	    
	    Arrays.setAll(lineX, i -> i * dis);
	    
	    String XYLineFile = dir + "Line100.dat";
	    Tecplot.lineXYTime(XYLineFile, lineX, lineY, 0.0);
	    
	    String heatFile = dir + "Heat100.dat";
	    mesh2D.toTecplot(heatFile, 0.0, x);
	    
	    String deformFile =  dir + "Deform100.dat";
	    mesh2D.toTecplotDeform(deformFile, 0.0, disp, 0.0);
	    
	    String expansionFile = dir + "expansion100.dat";
	    String temp00File = dir + "temp00_100.dat";
	    FileWriter fw1 = null;
	    FileWriter fw2 = null;
	    try {
			fw1 = new FileWriter(expansionFile);
			fw2 = new FileWriter(temp00File);
			fw1.write("0.0\t0.0\n");
			fw2.write("0.0\t0.0\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	    
	    double time = 0.0;
	    while(time < endTime) {
	    	time += dt;
	        System.out.println("time: " + time);
	        fs.applyBC_MBN(RHS, Direct.All, 1.0, 3);
	        solver.PCGSSOR(RHS, x, 1.5, 1);
	        RMatrix.mul(x, RHS);
	        
	        mesh2D.toTecplot(heatFile, time, x);
	        Arrays.setAll(lineY, i -> fs.getValue(x, new double[] {0.25, lineX[i]}));
	        Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);
	        
	        Arrays.fill(heatLoad, 0.0);
	        vecX[0] = x;
	        fs2.assembleHeatLoad(heatLoadFunc, null, vecX, BVPType.PLANESTRAIN, heatLoad);
	        //System.out.println(MVO.toString(heatLoad));
	        fs2.applyBCZero(heatLoad, Direct.X, 2, 4);
	        fs2.applyBCZero(heatLoad, Direct.Y, 1);
	        solver2.PCGSSOR(heatLoad, d, 1.5, 1);
	        fs2.extractComponent(d, disp);
	        mesh2D.toTecplotDeform(deformFile, time, disp, (1-Nu)/(alpha * (1 + Nu)));
	        
	        try {
	        	fw1.write(time + "\t " + fs.getValue(disp[1], new double[] {0.25, 1.0}) * (1-Nu)/(alpha * (1 + Nu)) + "\n");
				fw2.write(time + "\t " + fs.getValue(x, new double[] {0.0, 0.0}) + "\n");
			} catch (IOException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}
	    }
	    try {
			fw1.close();
			fw2.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	    

	}

}
